Глава III
3.3. Производные кривых B-сплайнов

Пусть `C^{(k)}(u)` обозначает `k`-ю производную `C(u)`. Если u фиксирован, мы можем получить `C^((k))(u)`, вычисляя `k`-ю производную базисных функций (смотри формулы [2.7], [2.9], и [2.10] и Алгоритм A2.3). в частности

`C^((k))(u)=sum_{i=0}^nN_{i,p}^((k) )P_i`(3.3)

Рассмотрим пример из раздела 2.5, при `p=2`, `U={0,0,0,1,2,3,4,4,5,5,5}`, и `u=`1⁄2. Из уравнения (2.7) мы имеем

`N_{2,2}^'(5/2)=0-2/(3-1)1/2=-1/2`

`N_{3,2}^'(5/2)=2/(3-1)1/2-2/(4-1)1/2=0`

`N_{4,2}^'(5/2)=2/(4-1)1/2-0=1/2`

Отсюда следует, что

`C^'(5/2)=-1/2P_2+1/2P_4`

Алгоритм для вычисления точки на кривой B-сплайна и все производные до и вклю­чая `d`ую по фиксированному значению u следующим образом. Мы позволяем `d>p`, хотя эти производные 0 в этом случае (для нерациональных кривых); эти прои­зводные являются необходимыми для рациональных кривых. Вход для алгоритма `u`, `d`, и кривая B-сплайна, определенный (в течение оставшейся части этой книги) с помощью

`n`∶ число контрольных точек `n+1`;

`p`∶ степень кривой;

`U`∶ узлы;

`P`∶ контрольные точки;

Выход это массив CK[], где CK[k] это `k`-ая производная, `0≤k≤d`. Мы используем Алгоритмы А2.1 и А2.3. Локальный массив, nders[][], используется для хранения производных базисных функций.

Алгоритм А3.2

CurveDerivsAlg1(n,p,U,P,u,d,CK)
{ /* Вычислить  производные  кривой */
  /* Вход: n,p,U,P,u,d  */
  /* Выход: CK  */
du = min(d,p);
for (k=p+1; k<=d; k++) CK[k] = 0.0;
span = FindSpan(n,p,u,U);
DersBasisFuns(span,u,p,du,U,nders);
for (k=0; k<=du; k++)
   {
   CK[k] = 0.0;
   for (j=0; j<=p; j++)
      CK[k] = CK[k] + nders[k][j]*P[span-p+j];
   }
}
    

Теперь вместо фиксированного `u`, мы хотим официально дифференцировать `p`-ую степень кривой B-сплайна,

`C(u)=sum_{i=0}^nN_{i,p}(u)P_i`

определяется по вектору узла

`U="{"ubrace(0,…,0)_(p+1),u_{p+1},…,u_{m-p-1},ubrace(1,…,1)_(p+1)"}"`

Из уравнений (3.3) и (2.7) мы получаем

`C^'(u)=sum_{i=0}^nN_{i,p}^'P_i` `=sum_{i=0}^n(p/{u_{i+p}-u_i}N_{i,p-1}(u)-p/{u_{i+p+1}-u_{i+1}}N_{i+1,p-1}(u))P_i` `=(psum_{i=-1}^{n-1}N_{i+1,p-1}(u)P_{i+1}/{u_{i+p+1}-u_{i+1}})` `-(psum_{i=0}^nN_{i+1,p-1}(u)P_i/{u_{i+p+1}-u_{i+1}})` `=p{N_{0,p-1}(u)P_0}/{u_p-u_i}+psum_{i=-1}^{n-1}N_{i+1,p-1}(u)((P_{i+1}-P_i))/{u_{i+p+1}-u_{i+1}}` `-p{N_{n+1,p-1}(u)P_n}/{u_{n+p+1}-u_{n+1}}`

Первый и последний члены выражаются в 0⁄0, что 0 по определению. Таким образом

`C^'(u)=psum_{i=-1}^{n-1}N_{i+1,p-1}(u)((P_{i+1}-P_i))/{u_{i+p+1}-u_{i+1}}=sum_{i=-1}^{n-1}N_{i+1,p-1}(u)Q_i`

где   `Q_i=p{P_{i+1}-P_i}/{u_{i+p+1}-u_{i+1}}`(3.4)

Теперь пусть `U^'` будет узловым вектором полученным путем отбрасывания первого и последнего узлов из `U`, т.е.

`U^'="{"ubrace(0,…,0)_(p),u_{p+1},…,u_{m-p-1},ubrace(1,…,1)_(p)"}"`(3.5)

(`U^'` имеет `m-1` узлов). Тогда легко проверить, что функция `N_{i+1,p-1}(u)`, рассчитанная на `U`, равна `N_{i,p-1}(u)` вычисляется на `U^'`. Таким образом

`C^'(u)=sum_{i=0}^{n-1}N_{i,p-1}(u)Q_i`(3.6)

где `Q_i` определены по формуле (3.4) и `N_{i,p-1}(u)` рассчитываются на `U^'`. Таким образом, `C^'(u)` является `p-1`й-степени кривая B-сплайна.

Примеры

Пример3.1 Пусть `C(u)=sum_{i=0}^4N_{i,2}(u)P_i` квадратичная кривая, определенная на

`U={0,0,0,`2⁄5`,`3⁄5`,1,1,1}`

Тогда `U={0,0,`2⁄5`,`3⁄5`,1,1}` и `C^'(u)=sum_{i=0}^3N_{i,1}(u)Q_i` где

`Q_0={2(P_1-P_0)}/{2/5-0}=5(P_1-P_0)`

`Q_1={2(P_2-P_1)}/{3/5-0}=10/3(P_2-P_1)`

`Q_2={2(P_3-P_2)}/{1-2/5}=10/3(P_3-P_2)`

`Q_3={2(P_4-P_3)}/{1-3/5}=5(P_4-P_3)`

`C(u)` и `C'(u)` показаны на рисунках 3.15a и 3.15b соответственно.

Рисунок 3.15. (a) квадратичная кривая на `U={0,0,0,`2⁄5`,`3⁄5`,1,1,1}`; (b) производная кривой представляет собой кривую B-сплайна первой степени на `U={0,0,`2⁄5`,`3⁄5`,1,1}`.

Пример3.2 Пусть `C(u)=sum_{i=0}^6N_{i,3}(u)P_i` будет кубической кривой, определенной на

`U={0,0,0,0,`2⁄5`,`3⁄5`,`3⁄5`,1,1,1,1}`

Тогда `U={0,0,0,`2⁄5`,`3⁄5`,`3⁄5`,1,1,1}` и `C^'(u)=sum_{i=0}^5N_{i,2}(u)Q_i` где

`Q_0={3(P_1-P_0)}/{1/3-0}=9(P_1-P_0)`

`Q_1={3(P_2-P_1)}/{2/3-0}=9/2(P_2-P_1)`

`Q_2={3(P_3-P_2)}/{1-1/3}=9/2(P_3-P_2)`

`Q_3={3(P_4-P_3)}/{1-1/3}=9/2(P_4-P_3)`

`Q_4={3(P_5-P_4)}/{1-2/3}=9(P_5-P_4)`

`Q_5={3(P_6-P_5)}/{1-2/3}=9(P_6-P_5)`

`C(u)` и `C'(u)` показаны на рисунках 3.16a и 3.16b соответственно. Обратите внимание, что `C'(u)` - квадратичная кривая с острием в двойном узле `u=`3⁄5.

Рисунок 3.16. (a) кубическая кривая для `U={0,0,0,0,`2⁄5`,`3⁄5`,`3⁄5`,1,1,1,1}`; (b) кривая квадратичной производной на `U={0,0,0,`2⁄5`,`3⁄5`,`3⁄5`,1,1,1}`.

Пример3.3 Напоминая, что кривая Безье `p`-й степени является кривой B-сплайна на `U={0,...,0,1,...,1}` (без внутренних узлов), уравнение (3.4) сводится к `Q_i=p(P_{j+1}-P_i)` для `0<=i<=n-1`. Так как `n=p` и `N_{i,p-1}(u)=B_{i,n-1}(u)`, полиномы Бернштейна, уравнение (3.6) эквивалентно уравнению (1.9).

Первые производные на концах кривой B-сплайна даются

`C^'(0)=Q_0=p/u_{p+1}(P_1-P_0)`

`C^'(1)=Q_{n-1}=p/{1-u_{m-p-1}}(P_n-P_{n-1})`(3.7)

(смотри примеры Пример3.1 и Пример3.2, и рисунки 3.15 (a) и (b) и рисунки 3,16 (a) и (b)). Обратите внимание, что на рисунках 3.15b и 3.16b векторы производных и контрольные точки различно масштабированны для лучшей визуализации, через 1⁄2 и 1⁄3, соответственно.

Поскольку `C^'(u)` является кривой B-сплайна, мы применяем формулы с (3.4) до (3.6) рекурсивно, чтобы получить более высокие производные. Принимая `P_i^((0))=P_i`, мы пишем

`C(u)=C^((0))(u)=sum_{i=0}^nN_{i,p}(u)P_i^((0))`

Тогда   `C^((k))(u)=sum_{i=0}^{n-k}N_{i,p-k}(u)P_i^((k))`(3.8)

где   `P_i^((k))={(P_i,k=0),({p-k+1}/{u_{i+p+1}-u_{i+k}}(P_{i+1}^((k-1))-P_i^((k-1))),k>0):}`

и   `U^((k))="{"ubrace(0,…,0)_(p-k+1),u_{p+1},…,u_{m-p-1},ubrace(1,…,1)_(p-k+1)"}"`

Алгоритм A3.3 является не рекурсивной реализацией уравнения (3.8). Он вычи­сляет контрольные точки всех производных кривых до и включая `d`-ую производную (`d≤p`). На выходе PK[k][i] является `i`-ой управляющей точкой `k`-ой производной кри­вой, где `0<k<d` и `r_1≤i≤r_2-k`. Если `r_1=0` и `r_2=n`, вычисляются все контрольные точки.

Алгоритм А3.3

CurveDerivCpts(n,p,U,P,d,r1,r2,PK)
{ /* Вычислить  контрольные  точки  производных  кривой */
  /* Вход: n,p,U,P,d,r1,r2  */
  /* Выход: PK  */
r = r2 – r1;
for (i=0; i<=r; i++)
   PK[0][i] = P[r1+i];
for (k=1; k<=d; k++)
   {
   tmp = p – k + 1;
   for (i=0; i<=r-k; i++)
      PK[k][i]=tmp*(PK[k-1][i+1]- PK[k-1][i])/(U[r1+i+p+1] – U[r1+i+k]);
   }
}
    

Используя уравнение (3.8), мы вычисляем вторую производную в конечной точке, `u=0`, кривой B-сплайна (`p>1`).

`C^((2))(0)=P_0^((2))={p-2+1}/{u_{p+1}-u_2}(P_1^((1))-P_0^((1)))` `={p-1}/{u_{p+1}-u_2}[p/{u_{p+2}-u_2}(P_2^((0))-P_1^((0)) )-p/{u_{p+1}-u_1}(P_1^((0))-P_0^((0)))]`

Из `u_1=u_2=0` следует следущее

`C^((2))(0)={p(p-1)}/u_{p+1}[P_0/u_{p+1}-{(u_{p+1}+u_{p+2})P_1}/{u_{p+1}u_{p+2}}+P_2/u_{p+2}]`(3.9)

Аналогично,

`C^((2))(1)=p(p-1)/{1-u_{m-p-1}}xx` `[P_n/{1-u_{m-p-1}}-{(2-u_{m-p-1}-u_{m-p-2})P_{n-1}}/{(1-u_{m-p-1})(1-u_{m-p-2})}+P_{n-2}/{1-u_{m-p-2}}]` (3.10)

Обратите внимание, что для кривых Безье эти уравнения сводятся к соответствующим выражениям уравнения (1.10). Рисунок 3.17 показывает квадратичную кривую Рисунок 3.15a с векторами `C^((2))(0)` и `C^((2))(1)`. `C^((2))(u)` представляет собой кусочную кривую нулевой степени, то есть, он является постоянным (но разным) вектором на каждом из трех отрезков [0,2⁄5), [2⁄5,3⁄5), и [3⁄5,1].

Рисунок 3.17. Вторые производные в конечных точках кривой на рисунке 3.15a.

Мы заканчиваем этот раздел с другим алгоритмом вычисления точки на кривой B-сплайна и всех производных до и включая `d`-ую производную для фиксированного зна­чения `u` (сравните с Алгоритмом A3.2). Алгоритм основан на уравнении (3.8) и Алго­ритм A3.3. Мы предполагаем использование, AllBasisFuns, которая является простой модификацией BasisFuns (Алгоритм A2.2), чтобы вернуть все ненулевые базисные функции всех степеней от 0 до `p`. В частности, N[j][i] является значением базисной функции `i`-ой степени, `N_{span-i+j,i}(u)`, где `0≤i≤p` и `0≤j≤i`.

Алгоритм А3.4

CurveDerivsAlg2(n,p,U,P,u,d,CK)
{ /* Вычислить  производные  кривой */
  /* Вход: n,p,U,P,u,d  */
  /* Выход: CK  */
du = min(d,p);
for (k=p+1; k<=d; k++)  CK[k] = 0.0;
span = FindSpan(n,p,u,U);
AllBasisFuns(span,u,p,U,N);
CurveDerivCpts(n,p,U,P,du,span-p,span,PK);
for (k=0; k<=du; k++)
   {
   CK[k] = 0.0;
   for (j=0; j<=p-k; j++)
      CK[k] = CK[k] + N[j][p-k]*PK[k][j];
   }
}
    
Рисунок 3.18 показывает кубическую кривую с первой, второй и третьей производной, вычисленных для `u=`2⁄5. (Производные уменьшены на 2⁄5).

Рисунок 3.18. Кубическая кривая на `U={0,0,0,0,`1⁄4`,`3⁄4`,1,1,1,1}` с первой, второй и третьей производными, вычисленными при `u=`2⁄5.